// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                              avtPoincareIC.C                             //
// ************************************************************************* //

#include <avtPoincareIC.h>

#include <list>
#include <cmath>
#include <cassert>
#include <iostream>
#include <limits>
#include <ImproperUseException.h>
#include <DebugStream.h>
#include <avtVector.h>
#include <algorithm>

#include <vtkPlane.h>

// ****************************************************************************
//  Method: avtPoincareIC constructor
//
//  Programmer: Hank Childs
//  Creation:   October 4, 2010
//
// ****************************************************************************

avtPoincareIC::avtPoincareIC(
    int maxSteps_,
    bool doTime_,
    double maxTime_,
    unsigned int mask,
    const avtIVPSolver* model, 
    Direction dir,
    const double& t_start,
    const avtVector &p_start,
    const avtVector &v_start,
    int ID) :
    avtStateRecorderIntegralCurve(mask, model, dir, t_start, p_start, v_start, ID)
{
    numIntersections = 0;
    maxIntersections = 0;

    terminatedBecauseOfMaxIntersections = false;

    puncturePeriod          = 0;
    puncturePeriodTolerance = 0;

    numSteps = 0;
    maxSteps = maxSteps_;

    doTime = doTime_;
    maxTime = maxTime_;

    terminatedBecauseOfMaxSteps  = false;
}


// ****************************************************************************
//  Method: avtPoincareIC constructor
//
//  Programmer: Hank Childs
//  Creation:   October 4, 2010
//
// ****************************************************************************

avtPoincareIC::avtPoincareIC() : avtStateRecorderIntegralCurve()
{
    numIntersections = 0;
    maxIntersections = 0;

    puncturePeriod          = 0;
    puncturePeriodTolerance = 0;

    terminatedBecauseOfMaxIntersections = false;

    puncturePeriod          = 0;
    puncturePeriodTolerance = 0;

    numSteps = 0;
    maxSteps = 0;

    doTime = false;
    maxTime = 0;

    terminatedBecauseOfMaxSteps  = false;
}


// ****************************************************************************
//  Method: avtPoincareIC destructor
//
//  Programmer: Hank Childs
//  Creation:   October 4, 2010
//
// ****************************************************************************

avtPoincareIC::~avtPoincareIC()
{
}


// ****************************************************************************
//  Method: avtPoincareIC::CheckForTermination
//
//  Purpose:
//      Checks to see if we should terminate the integral curve.
//
//  Programmer: Hank Childs
//  Creation:   October 4, 2010
//
//  Modifications:
//
//    Hank Childs, Sun Dec  5 10:18:13 PST 2010
//    Add avtIVPField argument.
//
// ****************************************************************************

bool
avtPoincareIC::CheckForTermination(avtIVPStep& step, avtIVPField *)
{
    bool shouldTerminate = false;

    // Check for termination.
    if (numSteps && IntersectPlane( step.GetP0(), step.GetP1(),
                                    step.GetT0(), step.GetT1() ))
    {
        numIntersections++;

        if (numIntersections >= maxIntersections)
        {
            terminatedBecauseOfMaxIntersections = true;
            shouldTerminate = true;
        }
    }

    // When doing a double puncture plot many points are examined and
    // one may not get enough points for the analysis thus allow a cut
    // off so not to integrate forever.
    if( puncturePeriod )
    {
      if( doTime )
      {
        if( (direction == DIRECTION_FORWARD  && step.GetT1() >= maxTime) ||
            (direction == DIRECTION_BACKWARD && step.GetT1() <= maxTime) )
            shouldTerminate = true;
      }
    }      
     
    ++numSteps;

    // If max steps is zero ignore the test.
    if( !shouldTerminate && maxSteps && numSteps >= maxSteps )
    {
      terminatedBecauseOfMaxSteps = true;
      shouldTerminate = true;
    }

    return shouldTerminate;
}


// ****************************************************************************
//  Method: avtPoincareIC::SetIntersectionCriteria
//
//  Purpose:
//      Defines an object for integral curve intersection.
//
//  Programmer: Dave Pugmire
//  Creation:   August 10, 2009
//
//  Modifications:
//
//   Dave Pugmire, Tue Aug 18 08:47:40 EDT 2009
//   Store plane equation.
//
// ****************************************************************************

void
avtPoincareIC::SetIntersectionCriteria(vtkObject *obj, int maxInts)
{
    // Only plane supported for now.
    if (!obj->IsA("vtkPlane"))
        EXCEPTION1(ImproperUseException, "Can not SetIntersectionCriteria in avtPoincare, the Poincare plot only supports plane intersections.");


    intersectPlanePt   = avtVector(((vtkPlane *)obj)->GetOrigin());
    intersectPlaneNorm = avtVector(((vtkPlane *)obj)->GetNormal());

    intersectPlaneNorm.normalize();
    intersectPlaneEq[0] = intersectPlaneNorm.x;
    intersectPlaneEq[1] = intersectPlaneNorm.y;
    intersectPlaneEq[2] = intersectPlaneNorm.z;
    intersectPlaneEq[3] = intersectPlanePt.dot(intersectPlaneNorm);

    maxIntersections = maxInts;
}


// ****************************************************************************
//  Method: avtPoincareIC::SetPuncturePeriodCriteria
//
//  Purpose:
//      Defines an object for integral curve puncture period intersection.
//
//      This criteria is used for the double Poincare section
//
//  Programmer: Allen Sanderson
//  Creation:   April 10, 2014
//
//  Modifications:
//
// ****************************************************************************

void
avtPoincareIC::SetPuncturePeriodCriteria( double period, double tolerance )
{
    puncturePeriod          = period;
    puncturePeriodTolerance = tolerance;
}


// ****************************************************************************
//  Method: avtPoincareIC::IntersectPlane
//
//  Purpose:
//      Intersect integral curve with a plane.
//
//  Programmer: Dave Pugmire
//  Creation:   August 10, 2009
//
//  Modifications:
//
//   Dave Pugmire, Tue Aug 18 08:47:40 EDT 2009
//   Don't record intersection points, just count them.
//
//   Dave Pugmire, Tue Dec  1 11:50:18 EST 2009
//   Switch from avtVec to avtVector.
//
// ****************************************************************************

bool
avtPoincareIC::IntersectPlane(const avtVector &p0, const avtVector &p1,
                              const double    &t0, const double    &t1)
{
    bool foundIntersection = false;

    double distP0 = intersectPlaneEq[0] * p0.x +
                    intersectPlaneEq[1] * p0.y +
                    intersectPlaneEq[2] * p0.z -
                    intersectPlaneEq[3];

    double distP1 = intersectPlaneEq[0] * p1.x +
                    intersectPlaneEq[1] * p1.y +
                    intersectPlaneEq[2] * p1.z -
                    intersectPlaneEq[3];

#define SIGN(x) ((x) < 0.0 ? -1 : 1)

    // If either point on the plane, or points on opposite
    // sides of the plane, the line intersects.
    if( // distP0 == 0.0 || distP1 == 0.0 ||
        SIGN(distP0) != SIGN(distP1) )
    {
      avtVector dir(p1-p0);
      
      double dot = intersectPlaneEq[0] * dir.x +
                   intersectPlaneEq[1] * dir.y +
                   intersectPlaneEq[2] * dir.z;

      // If the segment is in the same direction as the plane then
      // count the intersection.
      if( dot > 0.0 )
      {
        // Double Poincare puncture test.
        if( puncturePeriod > 0 )
        {
          avtVector w = p1 - intersectPlanePt;
          
          // Find where it intersects the plane.
          double t = -(intersectPlaneEq[0] * w.x +
                       intersectPlaneEq[1] * w.y +
                       intersectPlaneEq[2] * w.z -
                       intersectPlaneEq[3]) / dot;
          
          avtVector point = avtVector(p1 + dir * t);

          double time = t0 + (t1-t0) * (point - p0).length() / (p1-p0).length();

          // Get the number of periods traversed.
          double nPeriods = time / puncturePeriod;

          // Get the factional part - should be close
          // to zero for an even period.
          double intPart, fracPart = modf(nPeriods, &intPart);
          
          if( fracPart < puncturePeriodTolerance ||
              1.0-puncturePeriodTolerance < fracPart )
            foundIntersection = true;
        }
        else
            foundIntersection = true;

      }
    }
    
    return foundIntersection;
}


// ****************************************************************************
//  Method: avtPoincareIC::Serialize
//
//  Purpose:
//      Serializes a integral curve so it can be sent to another processor.
//
//  Programmer: Hank Childs
//  Creation:   October 4, 2010
//
//   David Camp, Wed Mar  7 10:43:07 PST 2012
//   Added a Serialize flag to the arguments. This is to support the restore
//   ICs code.
//
// ****************************************************************************

void
avtPoincareIC::Serialize(MemStream::Mode mode, MemStream &buff, 
                         avtIVPSolver *solver, SerializeFlags serializeFlags)
{
    // Have the base class serialize its part
    avtStateRecorderIntegralCurve::Serialize(mode, buff, solver, serializeFlags);

    buff.io(mode, intersectPlanePt);
    buff.io(mode, intersectPlaneNorm);

    buff.io( mode, intersectPlaneEq[0]);
    buff.io( mode, intersectPlaneEq[1]);
    buff.io( mode, intersectPlaneEq[2]);
    buff.io( mode, intersectPlaneEq[3]);

    buff.io(mode, numSteps);
    buff.io(mode, maxSteps);
    buff.io(mode, doTime);
    buff.io(mode, maxTime);
    buff.io(mode, terminatedBecauseOfMaxSteps);
    buff.io(mode, numIntersections);
    buff.io(mode, maxIntersections);
    buff.io(mode, terminatedBecauseOfMaxIntersections);
}

// ****************************************************************************
//  Method: avtPoincareIC::MergeIntegralCurve
//
//  Purpose:
//      Merge a values from one curve into another
//
//  Programmer: Allen Sanderson
//  Creation:   August 4, 2015
//
// ****************************************************************************

void
avtPoincareIC::MergeIntegralCurve(avtIntegralCurve *ic)
{
    avtPoincareIC *pic = (avtPoincareIC *) ic;

    numSteps = pic->numSteps;
    terminatedBecauseOfMaxSteps = pic->terminatedBecauseOfMaxSteps;

    numIntersections = pic->numIntersections;
    terminatedBecauseOfMaxIntersections =
    pic->terminatedBecauseOfMaxIntersections;
}
